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Standard techniques used to model chemically-reacting 
flows require an artificial viscosity for stability in the presence 
of strong shocks. The resulting shock is smeared over at least 
three computational cells, so that the thickness of the shock is 
dictated by the structure of the overall mesh and not the shock 
physics. A gas passing through a strong shock is thrown into 
a nonequilibrium state and subsequently relaxes down over 
some finite distance to an equilibrium end state. The artificial 
smearing of the shock envelops this relaxation zone which 
causes file chemical kinetics of the flow to be altered. This 
paper presents a method which can investigate these issues by 
following the chemical kinetics and flow kinematics of a gas 
passing through a fully resolved shock wave at hypersonic 
Mach numbers. A nonequilibrium chemistry model for air is 
incorporated into a spectral multi-domain Navier-Stokes solu- 
tion method. Since no artificial viscosity is needed for stabil- 
ity of the multi-domain technique, the precise effect of this 
artifice on the chemical kinetics and relevant flow features can 
be determined. 
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1. INTRODUCTION 

Since the inception of the design of the Space Shuttle more than 15 
years ago, there has been considerable progress in prediction methodolo- 
gies. We are beginning tp see increasing sophistication in chemically 
reacting flow methods which involve nonequilibrium and two-temperature 
chemical kinetics [1,2]. To successfully couple a nonequilibrium chemistry 
model with the flow kinematics, we must first fully understand the impact 
the numerical algorithm has on the problem at hand. Such an investigation 
has not been researched to date; this work addresses fundamental questions 
relevant to such understanding. 

The standard technique used to model chemically reacting or non- 
reacting flows are finite-difference or finite- volume methods. One common 
feature shared by both techniques is the artificial viscosity needed in the 
presence of strong shocks. The effect that the artificial viscosity has is to 
smear the shock over at least three computational cells. The size of these 
cells is not directed by the shock physics but by the structure of the overall 
mesh. 

When air passes through a strong shock wave, there follows a conver- 
sion of energy of directed translational motion into other forms of energy 
such as translation, rotation, vibration, and dissociation. Since the tempera- 
ture of a gas depends on the energy of translation, when energy is 
transferred into this form, the temperature of the gas rises. If, at a later 
time, a portion of the energy of translation is converted into other modes, 
the temperature falls. Because of differences in the magnitude of the time 
required for the different processes to absorb their equilibrium amount of 
energy, considering them to occur separately and consecutively results in a 
good approximation [5]. The variation of temperature from its free-stream 
value to its value following establishment of thermal equilibrium can, 
therefore, be represented as in the schematic given in Fig. 1 [3]. Station 
(1) represents conditions in the free stream ahead of the shock wave where 
the temperature is low enough that vibration of the molecules is not appre- 
ciably excited. Rotation of the molecules is, however, completely excited. 
Station (2) represents conditions in the free stream ahead of the shock 
wave where the temperature is low enough that vibration of the molecules 
is not appreciably excited. Rotation of the molecules is, however, com- 
pletely excited. Station (2) represents conditions behind but very close to 
the shock wave. The thermal jump across the shock wave throws the 
translational and rotational degrees of freedom temporarily out of equili- 
brium, but only a few molecular collisions are required to establish equili- 
brium not only of the translational but also of the rotational degrees of 
freedom. Adjustment of the vibrational degrees of freedom occurs between 
stations (2) and (3). This process requires a number of collisions (or a 
time or distance) that are large compared with that required for adjustment 
of translation and rotation, but are small compared with that required for 
establishment of dissociation equilibrium. At station (3), vibration is fully 
established but dissociation is only beginning. At station (4), sufficient 
time has elapsed for the dissociation of oxygen to reach equilibrium; and, 
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at station (5), the dissociation of nitrogen is in equilibrium and the tem- 
perature has reached its equilibrium value. These relaxation pathways are 
transition zones which initiate the chemical kinetics which follow. The 
pertinent questions are: Does the artificial smearing of the shock envelop 
the transition region; and if so, is the distribution of energies and chemical 
composition significantly altered so as to affect the chemical kinetics of the 
flow. 
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Fig. 1. Schematic of relaxation zones for nonequilibrium air. 


Fig. 2 is a plot of data accumulated from shock-tube experiments 
[4,5, 6,7]. Shock thickness and relaxation path lengths are plotted against 
Mach number in air. The data corresponds to shock thicknesses and relax- 
ation path lengths which exist at an altitude of 120,000 feet As can be 
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seen, die shock thickness at this altitude is within two or three orders of 
magnitude of these relaxation lengths. A typical calculation easily smears 
a shock over several orders of magnitude which is sufficient to obliterate 
this important region. Though it is possible with standard finite-difference 
techniques to resolve die shock sufficiently so that vibrational relaxation 
can be modeled, a calculation for more than just a very small region will 
be computationally to expensive. To make such a calculation practical for 
full-scale application the needed resolution must be attainable on a tract- 
able grid. At the very least, present day computations of chemically- 
reacting flows should be paralleled by fundamental studies of die effects 
factors inherent in the algorithm have on the relevant physics. 

Relaxation length Shock thickness 

— -Logan (1957) POOP Camac 0964). Alsmeyer (1975) 

ffffH Wray, Freeman (1963) 



Fig. 2. Comparison of relaxation path lengths and shock thicknesses; 
altitude 'b 120,000 feet. 

A methodology well suited for carrying out such a study is die spec- 
tral multi-domain technique [7]. One of the strongest properties of spec- 
tral methods is their exponential rate of convergence and high-order accu- 
racy on relatively coarse meshes [8]. Extensions of their application to 
include capturing regions of high gradient and modeling complex 
geometries has recendy been attainable due to the spectral multi-domain 
technique developed by Macaraeg and Streett [9]. The multi-domain 
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technique interfaces regions with independent discretizations, while main- 
taining spectral accuracy and an exponential rate of error decay. Adjoin- 
ing regions are interfaced by enforcing a global flux balance which 
preserves high-order continuity of the solution. This interface technique 
maintains spectral accuracy, even when mappings and/or domain sizes are 
radically different across the interface. Such advantages are needed when 
computing chemically reacting flows, which can display widely disparate 
temporal and spatial scales due to fast chemical reaction as opposed to the 
macroscopic scales associated with the hydrodynamics. 

2. MULTI-DOMAIN SPECTRAL TECHNIQUE 

A number of other spectral domain decomposition techniques have 
appeared in the literature. For example, the spectral element method which 
applies finite element methodology using Galerkin spectral discretization 
in the variational formulation within elements is a popular technique 
[10,11]. One drawback with this technique is that it utilizes a split 
Galeridn-collocation discretization which restricts its application to 
convection-diffusion problems for incompressible flows. The spectral ele- 
ment method in practice is used in a manner similar to classical finite ele- 
ment techniques: low-order internal discretization using many elements 
with no internal stretchings to improve resolution. In addition each ele- 
ment contains the same number of collocation points. Other domain 
decomposition techniques involve explicit enforcement of Cl continuity 
across the interface [12,13]. It is not clear how well these techniques per- 
form for strongly convection-dominated problems; the second author’s 
experience with such techniques [14] has shown them to be not entirely 
satisfactory. 

The multi-domain spectral technique overcomes the above restrictions 
by splitting the domain into regions with independent collocation discreti- 
zations, each of which preserve the advantages of spectral collocation, and 
allow the ratio of the mesh spacings between regions to be several orders 
of magnitude higher than allowable in a single domain. The partitioning 
results in three classes of points: the interior points where collocation of 
the equation is applied, boundary points where the physical boundary con- 
ditions are applied, and interface points which satisfy a global flux balance 
across adjoining sub-domain and preserve high-order continuity of the solu- 
tion. The global flux balance is imposed in terms of a spectral integral of 
the discrete equations across adjoining domains. In the multi-dimensional 
case this integration is performed along the coordinate line perpendicular to 
the interface. If adjoining domains have coordinate lines perpendicular to 
the interface which do not align, the solution from the domain with the 
lesser number of points are spectrally interpolated from the fine grid onto 
the coarse grid. To close the equation set, continuity of the solution and 
the flux is imposed at the interface. The total system of interface relations 
has a block structure- each block spanning a single interface between a 
pair of adjoining domains. Since the discretizations interior to each domain 
are solved uncoupled from each other and since the interface relation has a 
block structure, the solution scheme can be adapted to the particular 
requirement for each domain. Iterative and noniterative schemes have been 
developed to solve the combined interior/interface system. Details of the 
algorithm as well as comparisons between the present method, spectral ele- 
ment, and single domain collocation are given in ref. 7. 
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3. NUMERICAL MODEL 

The above technique will model the chemical kinetics and flow 
kinematics of a nonionized air mixture (02, N2, NO, O, and N) passing 
through a fully resolved shock wave, thus alleviating the need for artificial 
viscosity. The governing equations are the quasi-one dimensional Navier- 
Stokes equations [15], and species conservation equation [1]. The quasi- 
one dimensional form is used to provide an artifice for controlling the 
shock location in the physical space for this otherwise indeterminate prob- 
lem. The equations are nondimensionalized by dividing the state and tran- 
sport parameters by their dimensional free stream values. 

The viscosity of each of the individual species are calculated from a 
curve fit relation [16]. Similarly, curve fits are used to obtain specific 
heats, internal energies, and enthalpies [17,18]. The thermal conductivity 
of each specie is calculated from the Euken semi-empirical formula using 
the specie’s viscosity and specific heat. Appropriate mixture rules are next 
used to obtain the above transport properties of the mixture [19]. Experi- 
mental values of bulk viscosities, as obtained from acoustical inter- 
ferometry and related experiments, are taken from Truesdell [20]. 

In the present work, the diffusion model is limited to binary diffusion 
with the binary diffusion coefficients specified by the Lewis number. The 
value of the Lewis number used is 1.4. 

The temperature range under study will not exceed 8000 degrees kel- 
vin, for conditions at an altitude of approximately 190,000 feet. Therefore, 
ionization reactions, which occur at roughly 9000 degrees kelvin, are not 
included. The chemical reactions utilized for the nonionized air mixture 
are impact dissociation and exchange reactions. The seventeen reactions 
included in the present study can be found in ref. 1, which also lists ioniza- 
tion reactions. The constants needed to evaluate reaction rates are given in 
ref. 2. 

Initial conditions are obtained from a spectral multi-domain code for 
solution of the Navier-Stokes equations with equilibrium chemistry, written 
for the above problem. These governing equations may be found in ref. 
15. Transport properties are obtained in the manner previously discussed. 
The routines of ref. 17 generalized for air are used to obtain equilibrium 
concentrations. 

4. NUMERICAL ALGORITHM 

The multi-domain discretization involves three independent sub- 
domains, with the shock located in the center subdomain. Shock jump 
conditions are obtained by an iterative procedure to solve the Rankine- 
Hugoniot relations for real air. 

A direct inversion of the coupled system is utilized to obtain a fully 
implicit method. The conserved variables are written in delta form, and a 
pseudo time iteration using backwards Euler is utilized to obtain the steady 
state solution. Because of the large rank and ill-condition of the Jacobian 
matrix, iterative improvement of the Gaussian elimination solution was 
found to be required. 
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5. METHOD VERIFICATION 

Validity of the multi-domain Navier-Stokes algorithm is demonstrated 
by comparison with experiment. A low-density wind tunnel study of 
shock-wave structure and relaxation phenomena in gases was conducted by 
Sherman [21]. The experiment measured shock wave profiles recorded in 
terms of the variation in the equilibrium temperature of a small diameter 
wire oriented parallel to the plane of the shock, as the wire is moved 
through the shock zone. Free stream Mach number is 1.98. For this test 
case, a Navier-Stokes spectral multi-domain calculation is performed for a 
perfect gas with temperature dependent properties and a nonzero bulk 
viscosity corresponding to air [20]. A comparison with experimental tem- 
peratures normalized by the free stream temperature versus normalized dis- 
tance is given in Fig. 3. The experimental data points are represented by 
the open symbols. The numerical results fall within a symbol width of the 
data. The multi-domain technique utilized three domains. The center 
domain, located between x = -.15, and x = 0.3, contains 21 points; the 
outer domains contain 11 points each. The computational domain spans -1 
to 1. The unit Reynolds number of the flow is 80. A calculation for a unit 
Reynolds number of 1000 is given in Fig. 4, showing the ability of the 
method to accurately resolve strong gradients without numerical oscilla- 
tions. The plot is of Mach number versus normalized distance. Three 
domains are again used; the center domain contains 17 points, the outer 
domains contain 11 points each. The interfaces are located at -.15 and -.1. 



Fig. 3. Comparison of multi-domain Navier-Stokes calculation mrii 
experimentally obtained temperatures, ' M 00 = 1.98, Re — 80. 
Discretization 11121111, interfaces at -.15 and 3. 



Fig. 4. Computed solution for Re = 1000; discretization 11117111, 
interface at -.15 and -.1. 
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6. RESULTS 

The case to be discussed with respect to effects of nonequilibrium 
assumptions on the flow/chemistry is for Moo = 11.0, Too= 350K, and poo = 
6 x 10' 8 g/cm 3 . These conditions invoke primarily O 2 dissociations, with 
N 2 dissociations just beginning. Temperatures are not yet high enough for 
ionization to occur, so vibrational energy modes remain unexcited. 

Typical discretizations used in this study were 15, 27 and 33 points in 
the upstream, middle and downstream domains, respectively. The 
backward-Euler implicit time-stepping algorithm typically required less 
than 2000 iterations to converge from an equilibrium starting solution, with 
at least an eight order of magnitude reduction in maximum residual. 

It is interesting to note the effect of downstream boundary conditions 
of the chemical species on the numerical results. Two boundary conditions 
are compared in Fig. 5. The first involves imposing the equilibrium end 
state concentrations corresponding to the post shock temperature and den- 
sity values from the equilibrium jump conditions. The second involves a 
zero gradient boundary condition on species’ concentrations. The plot is of 
the log of the concentration in the relaxation zone versus normalized dis- 
tance behind the shock. It is evident that the extrapolation condition leads 
to a smooth profile of concentrations down to a desired end state, whereas 
an imposed equilibrium en d state causes oscillations in the numerics and a 
resulting discontinuity between the imposed equilibrium end state and the 
pathway chosen by the chemical kinetics calculation. 

" • " 10,1 - Otriehlet 

[NO] - Dlrlchlet 

[N] - Dlrlchlet 

— [oj - extrapolation 

[NO] - extrapolation 

— — [N] - extrapolation 
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To test the effect of the extrapolation boundary conditions on the 
chemical kinetics, a comparison is made between a computation for which 
the post-shock physical domain is long enough to allow the concentrations 
to relax to their final values, and one in which the domain is severly trun- 
cated. The results for [O 7 ] and [O] are shown in Fig. 6 , with the 
truncated-domain ( [- 1 , 1 ] ) solution represented by symbols, and the 
extended-domain ( [-1,120] ) solution shown as lines. As can be seen, the 
profiles are identical. 

o [0j] on [-1.1] 

0 [0] on [-1.1] 

[0 2 ] on [-1,120] 



Fig. 6. Comparison of [O 2 ] and [O] in the relaxation zone 
spanning [-1,1] and [-1,120]. 


The study of the effect of artificial viscosity on the flow physics was 
carried out by adding the equivalent of second-order artificial viscosity to 
the momentum, energy, and species’ concentration equations. The amount 
of artificial viscosity introduced was such that the shock was spread out to 
a thickness about three orders of magnitude wider than the fully-resolved 
no artificial viscosity solution. This width was chosen to represent the grid 
spacing of a typical shock-capturing computation on a full-scale 
configuration, such as that mentioned in the Introduction. Figures 7a. and 
7b. show the Mach number and temperature profiles for the resolved-shock 
and smeared-shock cases, respectively; note that the entire physical domain 
is shown in Fig. 7b, whereas only the near-shock region on a greatly 
expanded scale is plotted in Fig. 7a. For the resolved-shock case, end- 
points are at -1 and 200, and interface points are at -.3 and .1. The inter- 
face locations for the smeared shock case are at 65 and 100 , with endpoints 
at 0 and 270. The most important feature to note in comparing these 
profiles is the 20 % reduction in the temperature overshoot as a result of 
the artificial viscosity. We conjecture that this in not the result of simple 
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smearing, but is produced by the change in the relative magnitude of the 
spatial scales of the flow kinematics and chemical kinetics. The tempera- 
ture overshoot is produced by the conversion of kinetic energy to thermal 
energy, on a scale too small for the chemistry to respond. Thus the near- 
shock chemistry is essentially frozen, and the ratio of enthalpy to internal 
energy remains near its freestream value. The chemical reactions then 
begin to respond on their own scale, as dictated by the reaction rate con- 
stants used in the model. A fraction of this released thermal energy is 
absorbed by the reactions, lowering the temperature. When the shock is 
smeared, the release of thermal energy occurs on a scale nearer that on 
which the reactions operate, and therefore some of this energy is absorbed 
by the reactions while the release is still taking place. This reduces the 
overshoot in temperature. 

No Artificial Viscosity 




Fig. 7. Comparison of temperature and Mach number profiles for calculations 
a) without artificial viscosity, and b) with artificial viscosity. 
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Simple smearing of the shock from the artificial viscosity cannot 
account for this effect as can be realized by considering an equilibrium cal- 
culation. For such a calculation, only one spatial scale exists, and is dic- 
tated by the Reynolds number. Thus a change in scale or a change in Re 
are equivalent. Introduction of second-order artificial viscosity may, of 
course, be looked at as lowering the effective Re of a calculation. Thus 
any change in the shock profile by the introduction of artificial viscosity 
which cannot be accounted for by a simple scale change must be attributed 
to an alteration of the interplay between flow kinematics and chemical 
kinetic scales. This effect can be expected to strengthen with increasing 
Mach number. 

For the Mach 1 1 case considered here, although the path along which 
the chemistry relaxes is significantly altered in the near-shock region, the 
chemical end states from the computations with and without artificial 
viscosity are within 3-4% of each other. This can be seen in Figs. 8a, b 
and 9a, b, which show the profiles of [O 2 ] and [O] (Figs. 8a, b) and [N] and 
[NO] (Figs. 9a, b) in the relaxation zone. This is not to say, however, that 
this situation will always occur especially for higher Mach numbers where 
vibrational excitation (ionization) occurs. The reduction in the temperature 
overshoot could result in a large enough change in the relaxation path that 
the end state could be affected. For example, ionization reactions will not 
occur if artificial viscosity decreases this overshoot so that temperatures are 
not sufficiently high to initiate its onset 

7. DISCUSSION 

Although smooth, accurate and consistent solutions have been 
obtained for the species’ concentrations in the post-shock region, which is 
the region of interest here, we feel the solutions are not entirely satisfactory 
in the upstreams region and in the shock. The difficulty is that the range 
over which several of the species’ concentrations vary is too large to be 
accurately represented, given die Cray single-precision computations 
employed. Although iterative improvement is used in the Jacobian inver- 
sion of the time-stepping scheme to reduce the effects of finite word 
length, the usable dynamic range of the solution was limited to about six 
orders of magnitude. In the solution for [O], for example, the end state 
concentration is about 10' 2 ; the [O] solution in the free stream, however, 
oscillates at about ±_ 10‘ 8 , where the equilibrium concentration should be 
about 10* J \ This extremely small oscillation disappears in the shock 
region when reactions begin to take place, but it is still something of an 
annoyance. Fig. 3 of reference [1] also shows evidence of this difficulty in 
the near-shock region. We are investigating recasting species’ concentra- 
tion equations in terms of the logarithm of the concentrations, to better 
resolve the entire dynamic range of the solution. 
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